Random field Ising model on networks with inhomogeneous connections 



o 
o 

(N 
D 
OO 



Sang Hoon Lce/'Q Hawoong Jeong,^'Q and Jae Dong Noh^'0 

^Department of Physics, Korea Advanced Institute of Science and Technology, Daejeon 305-701, Republic of Korea 
^Department of Physics, University of Seoul, Seoul 130-743, Republic of Korea 

(Dated: February 6, 2008) 

We study a zero-temperature phase transition in the random field Ising model on scale-free net- 
works with the degree exponent 7. Using an analytic mean-field theory, we find that the spins are 
always in the ordered phase for 7 < 3. On the other hand, the spins undergo a phase transition 
from an ordered phase to a disordered phase as the dispersion of the random fields increases for 
7 > 3. The phase transition may be either continuous or discontinuous depending on the shape 
of the random field distribution. We derive the condition for the nature of the phase transition. 
Numerical simulations are performed to confirm the results. 
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I. INTRODUCTION 

Recently physicists have recognized that the underly- 
ing graphs or networks for interacting systems have an 
intriguing structure. Such complex networks are distinct 
from the periodic lattices in Euclidean space in many 
aspects. The structural property of the complex net- 
works has been studied intensively 0, 13 ■ Besides the 
network topology itself, traditional topics of statistical 
physics of complex networks have been investigated as 
well UlllllllllillSllllIillllll. statistical physi- 
cal systems of networks are attractive since they display 
theoretically interesting critical phenomena jlO, fill Il3 | . 
They are also attractive for a possible application to var- 
ious phenomena in social systems having a complex un- 
derlying network structure 0, 0, 0| . 

In contrast to the periodic lattices in Euclidean space, 
complex networks have an inhomogeneous structure. Re- 
cent studies reveal that the structural inhomogeneity 
plays an important role in critical phenomena on com- 
plex networks [1 Ellll El El Ei . However, in the 
study of the critical phenomena, systems with quenched 
disorder have received little attention with only a few 
exceptions 0, 0] . In this work we investigate the ef- 
fect of quenched disorder and structural inhomogeneity 
on the nature of a phase transition. For this purpose, we 
study a disorder-driven phase transition in the random 
field Ising model (RFIM) on scale-free (SF) networks. A 
SF network is characterized by the power-law degree dis- 
tribution Pd{k) ~ k~'^ with the degree exponent 7. The 
power-law distribution indicates that SF networks have 
an inhomogeneous structure and the degree exponent 7 
determines the strength of the inhomogeneity in struc- 
ture. 

The RFIM has attracted much attention in statisti- 
cal physics [Ullllllllllllllllllllilll. Being 



compared with the spin glass model where the quenched 
disorder is present in the interaction among spins |28| . 
the RFIM has a quenched random external magnetic 
field applied to each site. The quenched disorder leads 
to a phase transition from an ordered ferromagnetic 
phase to a disordered paramagnetic phase. In spite of 
the simpler structure of the RFIM than the spin glass 
model, there are still remaining questions and contro- 
versies, especially over the nature of the phase transi- 
tion [laEalSllilElil. 

We study the zero-temperature phase transition in the 
RFIM using an analytic mean-field theory. Our analy- 
sis shows that the shape of the random field distribution 
and the degree exponent 7 determine the nature of the 
disorder-driven phase transition. We also perform nu- 
merical simulations, which confirm the analytic results. 
This paper is organized as follows: In Sec. II, an analytic 
approach based on mean-field theory and its prediction 
of the phase transition nature are provided. Numerical 
simulations follow in Sec. HI, and Sec. IV summarize our 
work. 
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II. ANALYTIC MEAN-FIELD THEORY 



The Hamiltonian of the RFIM on a network is given 
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where Si = ±1 is the Ising spin variable of node i — 
1,2, . . . , N, J > is the ferromagnetic coupling strength 
between neighboring spins, and a^- is the adjacency ma- 
trix element of the network. The matrix element takes 
the value of 1 (0) if two nodes i and j are (not) connected 
via a link. The degree of a node i is given by ki = Uij . 
We will set J = 1 hereafter for notational simplicity. Here 
the external magnetic field hi is a quenched random vari- 
able, which is distributed identically and independently 
according to a distribution function p{h). We only con- 
sider a symmetric distribution - that is, p{h) = p(—h). 
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It is convenient to write 



(i). 



(2) 



where po{x) is a normalized [J po{x)dx ~ 1] function de- 
termining the shape of the distribution and A is a mea- 
sure of the disorder strength. 

The ferromagnetic coupHng J favors the ordered state 
with all spins up or down. The external magnetic field, 
however, tends to pin each spin to a random direction. 
The competition between them may lead to a phase tran- 
sition, which will be investigated at zero temperature us- 
ing mean-field theory. In the framework of the mean-field 
theory, each spin Si is assumed to be in equilibrium un- 
der the effective magnetic field hi = aijUij -\- hi where 
rUj = (sj) is the average local magnetization at node j. 
Hence, the average local magnetization should satisfy the 
coupled mean-field equation 



{si) — tanh 



-(3h, 



(3) 



where [3 = 1/ (ksT) with the Boltzmann constant ks and 
temperature T. 

Instead of solving Eq. ^ directly, we make a sim- 
plification by assuming that the local magnetization de- 
pends only on the degree and the magnetic field - that 
is, rrii = m{ki,hi). It could be valid when the network 
has no internal structure and all nodes with the same 
degree and the same random field are statistically equiv- 
alent Then m(fc, h) should satisfy 



m(fc, h) = tanh 



dh'p{h')ra{k' , h')kPd{k'\k) + (3h 

(4)" 

where Pd{k'\k) is the conditional probability that a neigh- 
borhood of a node with the degree k has the degree k' . 
The conditional probability Pd{k'\k) measures a correla- 
tion between degrees of adjacent nodes. Although many 
real-world networks display a nontrivial degree correla- 
tion |29j | , we focus our attention on uncorrelated networks 
in this work for analytic tractability. It will be interest- 
ing to study the effect of the degree correlation on critical 
phenomena, which we leave for future work. Without the 
correlation, the conditional probability is given by 



Pd(k'\k) = k'Pdik')/k, 

with the mean degree k 

Now we define the order parameter 
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kPdjk) 
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dh p{h)'m{k, h) 



(5) 
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as the weighted average of the local magnetization. Us- 
ing Eqs. and ©, one finds that the order parameter 
should satisfy 



kPdik) 



dk^ 



dh p{h) tanh(/3mfc + (3h) (7) 



in the continuum limit. We are interested in the zero- 
temperature limit where /? — > -l-oo. Using p{h) = p{—h) 
and p{h) = po{h/ A) / A, we finally obtain the self- 
consistency (SC) equation for the order parameter at zero 
temperature given by 



where 



m = /(to) 



G{x) 



dk'^^^G{km/A), 
k 



dx' pq{x'). 



(8) 



(9) 



The SC equation depends on the network inhomogene- 
ity through Pd{k), the shape of the random field distri- 
bution through G{x), and the disorder strength A. SF 
networks have the power-law degree distribution. We use 
the following explicit form for the degree distribution for 
further analysis: 



Pd{k) = ck- 



(10) 



for k > fcp. Here /cq is a cutoff and c = (7 — 1)A:J ^ is a 
normalization constant. The Ith moment of the degree, if 
it exists, will be denoted as kK As for the magnetic field 
distribution, we assume that ^0(2;) in Eq. |5J is analytic 
at X = 1^. Then, the function G{x) in Eq. © can be 
expanded as 
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G{x) = J2^' 

n=Q 
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where 61 = 2po(0), ^3 = pg(0)/3, and so on. It has the 
limiting behavior that G{x ^ 0) = and G{x 00) = 1. 

The phase transition nature is determined by the lead- 
ing behavior of /(m) near m = (see Fig. 0. For 
a bounded degree distribution - e.g., Poisson distribu- 
tion - one can insert Eq. Hll|l into Eq. ||SJ) and ex- 
pand /(m) into a series of m with odd-integer powers. 
However, with the power-law degree distribution, the 
function /(m) has a singular expansion. In that case, 
one needs to split the function G{x) into two parts as 
Gix) = Gr{x) + G,{x) where Gr{x) = E„<s 
and Gs{x) = G{x) — Gr{x). Here n is the largest integer 
among all satisfying 2n -f 1 < 7 — 2. We will show that 
Gr (Gs) contributes to f{m) a regular (singular) part 
consisting of integer (noninteger) powers of m. Since the 
leading behavior of f{m) depends on the value of n, we 
consider the following three cases separately. 



A. 7 > 5 case 

In this case, n > 1 and /(to) can be expanded as 
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FIG. 1: Schematic plots of f{m) (solid line), which has the 
infinite slope at m = (a), is convex at m = (b) and is 
concave at m = (c). The dashed line represents the graph 
of m, and the solid circle represents the solution for the SC 
equation m — f{m). 

where C2n+i = b2n+ik^"'^^ /k. Note that Gs{x) = 
as a; ^ and Gs{x) = 0{x^^-^) as x ^ oo. 
This property guarantees that the integral in the second 
term converges to a finite value. Hence we find that 

/,„,.2E2|E(^) + f«|!!(^)%c,K), (13) 

with 9 = min{7 — 2, 5}. 

The function f{m) has a finite slope at m = and 
changes its convexity depending on the shape of the ran- 
dom field distribution given hy pq(x). This property leads 
to the following conclusion: For Pq(0) > 0, f{m) is con- 
vex as shown in Fig. ^b). Then, the order parameter 
m jumps from a nonzero value to zero at a threshold 
value of A. That is to say, the system undergoes a first- 
order phase transition. For Pq(0) < 0, /(m) is concave 
as shown in Fig.^c). Therefore, the system undergoes a 
continuous phase transition at Ac = 2po(0)fc^//c and the 
order parameter scales as 

m^{A,-Af, (14) 

with the order parameter exponent 

/3 = l/2. (15) 

These results coincide with those for the mean-field 
model where all spins interact with all others | T^l20ll26| . 

B. 3 < 7 < 5 case 

In this range of 7, n = 1 and Gr{x) = 61a;. So the 
function /(m) is given by 

/M-^m+cQ^-7°° dxx^^^GAx). 

k VA/ \A/ Jrnko/A 

Note that Gs{x) = 0{x^) as x ^ and Gs{x) = 0{x) 
as x ^ 00 since Gs{x) = G{x) — bix. This property 
guarantees that the integral converges to a finite value. 
So we find that 

nm)='j^Q+cDQ'-\0im% (16) 



where the constant D is given by 

D = j dx x^-^{G{x) -2pq{Q)x). (17) 

The function f{m) has a finite slope at m = and 
changes its convexity depending on the sign of the con- 
stant D. This leads to the following conclusion: For pos- 
itive D, the system undergoes a first-order phase transi- 
tion. For negative the system undergoes a continuous 
phase transition at Ac = 2po{0)k^ /k and the order pa- 
rameter scales as in Eq. (|14|l with the critical exponent 

P=^- (18) 

7-3 

We want to stress that the transition nature is deter- 
mined by the whole shape of the random field distribution 
given by the function pq{x). For 7 > 5, it is determined 
by the sign of Po(0) which is related to the local shape 
of po{x) near a; = 0. On the contrary, it is the sign 
of the constant D that determines the transition nature 
for 3 < 7 < 5. Hence, one may have the continuous 
transition even with the magnetic field distribution with 
Pq{0) > and vice versa. 

One may have a negative D for a distribution p{h) 
which has a peak at h — and decreases monotonically 
as \h\ increases. In such a case, the Ising spins become 
disordered gradually as the disorder strength grows. On 
the other hand, one may have a positive D for a distribu- 
tion p{h) which has a peak at nonzero values oi h = ±ho 
and a deep valley at h = 0. In such a case, the random 
field breaks the order abruptly. 

C. 2 < 7 < 3 case 

In this range of 7, n = and Gs{x) = G{x). By chang- 
ing the integration variable A: to a; = mk/A in Eq. ||SJ|, 
and using Eq. (|10|l . we can write the integral as 

f^rn)^c(-) / dxx'-^G{x). (19) 

Jmko/A 

Note that G{x) vanishes (at most) linearly as x ^ and 
saturates to 1 for x ^ 1. These properties guarantee 
that the integral converges to a finite value in the limit 
TO ^ 0, which yields that 

/(m) = c'(to/A)^-2 + ©(m^), (20) 

with a constant c' — c dx x^~''G{x). The function 
/(to) has the infinite slope at m = and the SC equation 
TO = /(m) has a nonzero solution 

m ~ A-(^-2)/(3-7) (21) 

at all values of A [see Fig. ^a)]. Therefore, the sys- 
tem is "always magnetized" irrespective of the shape of 
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the random field distribution and the disorder strength. 
The SF network with 2 < 7 < 3, where the second mo- 
ment of degree diverges, is famous for its pecuhar behav- 
ior such as the absence of the percolation and epidemic 
threshold 0, 0, Isil Is^ . This "absence of magnetiza- 
tion threshold" is another example of such characteristic 
behavior. 

In summary, we have a general criterion for the nature 
of the zero-temperature phase transition of the RFIM 
with the symmetric random field distribution p{h) on SF 
networks with the degree distribution Pd{k) ~ k '' with- 
out the degree correlation. For 2 < 7 < 3, the system 
is always magnetized and there is no phase transition. 
For 7 > 3, the system displays a phase transition at a 
finite value of A. The transition may be either the first- 
order or continuous phase transition. The condition for 
the first-order transition is that D > [see Eq. IjlTfl ] or 
Pq{0) > for 3 < 7 < 5 or 7 > 5, respectively. In the 
opposite case the transition is the continuous one and 
the critical exponent for the order parameter is given by 
Eq. H15|) or (|18|) . respectively. Finally we add a remark 
that a logarithmic correction appears when 7 = 3 or 5. 



III. NUMERICAL SIMULATION 

We perform a numerical study of the RFIM with the 
Hamiltonian in Eq. ^ at zero temperature to confirm 
the analytic result. First we generate a SF network of 
N nodes and K = 2k N links with the degree exponent 
7 using the so-called static model [s^- The static model 
has no degree correlation except for the region where 
2 < 7 < 3 [i^l- Nevertheless, the "always magnetized" 
characteristic of the system from the mean-field analysis 
still holds for 2 < 7 < 3, as we will see. Random mag- 
netic fields are then assigned to each node according to 
a distribution function p{h) = po{h/ A) / A. The ground- 
state spin configuration {si} is found and the weighted 
order parameter 



m 



1 



N 



is calculated, where ki is the degree of node i. The or- 
der parameter is averaged over different samples to yield 
(to) . Note that each sample has a different realization of 
a network configuration and a different realization of ran- 
dom fields. The average over these samples corresponds 
to the order parameter defined in Eq. ©. The exact 
ground state of the RFIM can be found numerically by 
adopting the map ping of the RFIM onto the maximum 
flow problem (SJl- For details of the mapping and the 
numerical algorithm solving the problem, we refer reader 
to Ref. [33. 

As for the random field distribution p{h) = 
Po{h/ A)/ A, we use the two functions po{x) — p+{x) and 
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FIG. 2: (m) versus A with the magnetic field distribution 
p{h) = p+{h/A)/A and 7 = 2.5 (a), 4.0 (b), and 6.0 (c). The 
dashed line in (a) has a slope —1. 



Po (x) — p- {x) which are given by 



p-{x) 



TT (■KX\ 



(22) 
(23) 



in the interval —\<x<\ and zero outside the interval. 
These functions have the following properties: ?3+(0) > 
and £> > for p+(a;), and p" (0) < and Z? < for 
p^(x). So we can test the analytic result with these two 
distribution functions. 



A. Numerical result with p(/i) = p+{h/ A) / A. 

We present the numerical data for the sample averaged 
magnetization (to) in Fig. 12 They were obtained from 
the static model networks with 7 = 2.5, 4.0, and 6.0 of 
sizes N = 1000, . . . , 64000. 

At 7 = 2.5, the ferromagnetic order with nonzero to 
persists at high values of A. Moreover, the log-log plot 
in Fig. EJa) suggests that the magnetization decreases 
algebraically. This behavior is consistent with the an- 
alytic result TO - A"(''"2)/(3-7) Eq. (EU- Accord- 
ing to it, the decay exponent should be —1 at 7 = 2.5. 
There is a little discrepancy in the decay exponent. We 
attribute the apparent discrepancy to a finite-size effect 
since the decay exponent approaches the analytic result 
as N increases. However, we cannot exclude a possibility 
that it could be due to the negative degree correlation at 
7 = 2.5. 

Figures |2Ib) and[2Ic) show that the order parameter 
vanishes abruptly at a certain threshold of A, which is 
a characteristic of a first-order phase transition. In or- 
der to prove the first-order nature we study the order 
parameter histogram H{m) near the threshold. Numeri- 
cally the histogram is measured by the fraction of samples 
whose order parameter value lies between m and to-|-<5to, 
which is equal to H{m)5m. In Fig.|21 we present the his- 
togram H{m) obtained numerically on the SF networks 
of = 64000 nodes with Sm = 0.01. At small values 
of A the histogram is peaked at a nonzero value of m, 
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FIG. 3: The histogram on an arbitrary scale of the magne- 
tization at several values of A near the transition on the SF 
networks of A'' = 64000 nodes with 7 = 4.0 (left column) and 
6.0 (right column). 



while it is peaked at m = at high values of A. In the 
intermediate values of A, there appear two peaks in the 
histogram, which indicates the phase coexistence. The 
two-peak structure near the threshold confirms the first- 
order transition nature. 



B. Numerical result with p{h) =p_(/i/A)/A 

We present the numerical data obtained with the ran- 
dom field distribution p{h) = p^{h/A)/A in Fig. E| At 
7 = 2.5 [Fig. ^a)], the order parameter remains finite 
and decreases algebraically as A, which is consistent with 
Eq. At 7 = 4.0 and 6.0 [Figs. Efb) andEfc)], the 

order parameter shows a threshold behavior. Unlike the 
case with p{h) — p+(/i/A)/A, the order parameter ap- 
proaches zero smoothly as A increases. It indicates that 
the transition could be a continuous transition. 

In order to examine the transition nature, we measure 
the Binder parameter 



U = 1 



3(m2)2- 



(24) 



The Binder parameter is supposed to take a nontrivial 
value at a critical point with scale invariance. It takes 
a trivial value 2/3 and in an ordered phase and in a 
disordered phase, respectively, in the N 00 limit. A 
critical point will manifest itself as a crossing point in the 
plot of U versus A at different system sizes N. 

Figure |31 shows the Binder parameter for the three 
cases, each of which corresponds to2<7<3,3<7<5, 
and 7 > 5. In Fig. Eta), there is no crossing point and 
the value U = 2/3, corresponding to the ordered state, 
persists as the system size grows. This behavior clearly 
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FIG. 4: (m) versus A with p{h) = p_(/i/A)/A and 7 = 
2.5 (a), 4.0 (b), and 6.0 (c). The dashed line in (a) has a 
slope —1. 



shows that the system is always magnetized in the ther- 
modynamic limit. Figures |5lb) andlSJc) show that there 
appear the crossing points at Ac ~ 38.0 for 7 = 4.0 and 
Ac ~ 31.5 for 7 = 6 where the Binder parameter is scale 
invariant and the system is critical. It indicates that the 
transition is a continuous transition. 

Since the phase transition is a continuous one, we ex- 
pect that the order parameter satisfies the critical finite- 
size-scaling form 



(m) = iV-''/'''F((Ac - AY'n), 



(25) 



where /3 is the order parameter exponent and v' is the 
finite-size-scaling exponent. The scaling function F(x) 
has the limiting behavior that F(x — > 0) ~ const and 
F{x 00) x^^/"' so that 



(m) ~ (Ac - Af 
for > (Ac - A)-''' and 

(to) - N-f^/"' 



(26) 



(27) 



for A^ < (Ac - A)-'''. 

The finite-size-scaling form is used to obtain the criti- 
cal exponents /3 and v' . In Fig. El we present the scaling 
plot for the order parameter according to Eq. H25() with 
the exponent values that give the best data collapse. We 
estimate that /3 = 0.75 and i^' = 3.47 for 7 = 4.0 and 
/3 = 0.45 and ^ 2.81 for 7 = 6.0. 

Repeating the same analysis, we obtained the values 
of (3 at several values of 7 > 3. The numerical results 
are compared with the analytic mean-field results [see 
Eqs. ifT^ and (fT^ J in Fig.d For 7 > 5, the critical expo- 
nent P is indeed 1/2 from the simulation result. One finds 
that the values of /3 for 3 < 7 < 5 deviate slightly from 
the mean-field prediction. Although we suspect that this 
may be due to the singular dependence of /3 = 1/(7 — 3) 
near 7 = 3, we cannot exclude a possibility that it may be 
due to a limitation of the our mean-field approximation. 
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FIG. 5: U versus A for various system sizes with the degree 
exponent (a) 7 = 2.5, (b) 7 = 4, and (c) 7 = 6. 



IV. DISCUSSION AND CONCLUSIONS 

The main result of our work is that the RFIM on scale- 
free networks is always magnetized for 2 < 7 < 3 and it 
has a disorder-driven phase transition for 7 > 3 whose 
nature depends on the shape of the random field dis- 
tribution p{h) = po{h/A)/A. As for the nature of the 
transition, similar results have been known in regular 
lattices in high-dimensional Euclidean space 
the disorder-driven zero-temperature phase transition is 
first order or continuous for a convex (p"(0) > 0) or con- 
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FIG. 7: The critical exponent /3 as a function of the degree 
exponent 7. The solid line corresponds to the mean-field re- 
sult /3 = 1/(7 - 3) for 3 < 7 < 5 and /3 = 1/2 for 7 > 5. 
The numerical simulation is based on the finite-size-scaling 
method. 



cave (p"(0) < 0) random field distribution, respectively. 
Our result shows that the same criterion is valid for SF 
networks with 7 > 5. For 3 < 7 < 5, the criterion is 
replaced by positivity or negativity of the quantity D de- 
fined in Eq. H17|) . Roughly speaking, the distributions 
highly peaked at /i = give rise to the continuous 
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transition, while distributions highly peaked at nonzero 
h = zb/ig give rise to the first-order transition. One may 
have a distribution with Pg(0) > but with D < 0. For 
example, po(-l < a: < 1) = 3[a + (1 - a)x^]/[2{l + 2a)] 
with a = 3/4 is such a function. We checked numerically 
that it indeed leads to the continuous phase transition at 
7 = 4. 

In summary, we have investigated the RFIM on scale- 
free networks with inhomogeneous connections. The net- 
work topology, especially the degree exponent 7, is shown 
to affect the phase transition and the critical exponent. 
The shape of the random field distribution is also respon- 



sible for the nature of the phase transition. 
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